###



pacman::p_load(tidyverse, corrr, knitr, kableExtra, haven, fixest, dotwhisker, RColorBrewer, broom)

#load data on EA and prepare for analysis


ead <-  readRDS("Data/Regimedata/Measures_merged.rds")


names(ead)
##create lead values for different DVs
ead <- ead %>%
  group_by(Final_Code) %>% arrange(year) %>% mutate(F_RegType_lied_EA = dplyr::lead(RegType_lied_EA, 1),
                                               F_RegType_RoW_EA = dplyr::lead(RegType_RoW_EA, 1),
                                               F_Politytype_Anocracy = dplyr::lead(Politytype_Anocracy, 1),
                                               F_status_fh_PF = dplyr::lead(status_fh_PF, 1),
                                               F_AnckarRegtype_MP_Autocracy = dplyr::lead(AnckarRegtype_MP_Autocracy, 1),
                                               F_HTW_RegType_MP_Autocracy = dplyr::lead(HTW_RegType_MP_Autocracy, 1,),
                                               F_RegType_magaloni_EA = dplyr::lead(RegType_magaloni_EA, 1))










##create indicator of how many CW changes 

varlist <- c("RegType_lied_n", "RegType_RoW_n", "Politytype_n", "fh_status_n", "AnckarRegtype_n", "HTW_RegType_n", "RegType_magaloni_n")

ead$c <- 1



for (i in 1:length(varlist)) {
  
  variable <- varlist[i]
  var <- sym(variable)
  changevar <- sym(paste0("change_", variable))
  leadvar <- sym(paste0("F_", variable))
  periodvar <- sym(paste0(variable, "_period"))
  durationvar <- sym(paste0(variable, "_period", "_duration"))
  durationvar_ln <- sym(paste0("ln_", variable, "_period", "_duration"))
  
  ead <- ead%>%mutate(!!leadvar := dplyr::lead(!!var, 1))
  print("Lead Variable Created")
  
  
  ead <- ead%>%mutate(!!changevar := if_else(!!var != dplyr::lead(!!var, 1), 1, 0))
  print("Change Variable Created")
  
  ead <- ead %>%
    group_by(Final_Code) %>% arrange(year) %>%
    mutate(!!periodvar := cumsum(!!changevar))
  
  print("Period Variable Created")
  
  
  
  # create period duration (current)
  ead <- ead %>% arrange(year) %>%
    group_by(Final_Code, !!periodvar)%>%
    mutate(!!durationvar := cumsum(c),
           !!durationvar_ln := log(!!durationvar + 1))
  
  print("Duration Variables Created")
  
  ead <- ead%>%ungroup()
  
}              


names(ead)
 ead %>% select(Final_Code, year, RegType_lied_n, F_RegType_lied_n, change_RegType_lied_n, RegType_lied_n_period, RegType_lied_n_period_duration , ln_RegType_lied_n_period_duration)



# load Kim's data
df <- read_dta("Data/Replication/Kim/nam_data.dta")

names(df)



data <- left_join(df, ead %>%rename(ccode = Final_Code), by = c("year", "ccode"))

#save(data, file = "data/Merged_Measures_Kim.rda")

###Create Regional Averages

pacman::p_load(tidyverse, corrr, knitr, kableExtra, haven, fixest, dotwhisker, RColorBrewer, cshapes, maptools, rgeos)



data <- data%>%mutate(neigh_EA_lied = NA,
                      neigh_DEM_lied = NA,
                      neigh_EA_RoW = NA,
                      neigh_DEM_RoW = NA,
                      neigh_EA_Polity = NA,
                      neigh_DEM_Polity = NA,
                      neigh_EA_FH = NA,
                      neigh_DEM_FH = NA,
                      neigh_EA_AF = NA,
                      neigh_DEM_AF = NA,
                      neigh_EA_HTW = NA,
                      neigh_DEM_HTW = NA,
                      neigh_EA_AWD = NA,
                      neigh_DEM_AWD = NA)%>%ungroup()









yearlist <- 1960:2006
countrylist <- unique(data$ccode)




for(j in 1:length(yearlist)){
  
  
  current_year <- yearlist[j]
  when <- paste0("31/1/", current_year)
  
  
  print(paste("Loading distance matrix for", current_year))
  
  a <- distmatrix(
    date = as.Date(when,format='%d/%m/%Y'),
    type = "mindist",
    keep = 0.1,
    useGW = F,
    dependencies = FALSE
  )
  
  print(paste("Loaded distance matrix for", current_year))
  
  
  a <- as.data.frame(a)
  a$ccode <- rownames(a)
  
  
  
  
  
  
  
  for (i in 1:length(countrylist)){ 
    
    
    current_country <- countrylist[i]
    
    if(current_country %in% names(a)){
      
      print(paste("Beginning with:", current_country, current_year))
      
      
      countries <- data%>%distinct(ccode)
      
      
      
      b <- a%>%select(ccode, md = as.character(current_country))%>%filter(md <= 501 & ccode != all_of(current_country))
      
      c <- data%>%filter(ccode %in% unique(b$ccode) & year == current_year)%>%select(RegType_lied_EA, RegType_lied_DEM,
                                                                                     RegType_RoW_DEM, RegType_RoW_EA,
                                                                                     Politytype_Anocracy, Politytype_Democracy,
                                                                                     status_fh_F, status_fh_PF,
                                                                                     AnckarRegtype_MP_Autocracy,AnckarRegtype_Democracy,
                                                                                     HTW_RegType_MP_Autocracy, HTW_RegType_Democracy,
                                                                                     RegType_magaloni_EA, RegType_magaloni_DEM)
      
      
      
      
      EA_sh_lied <- mean(c$RegType_lied_EA, na.rm = T)
      DEM_sh_lied <- mean(c$RegType_lied_DEM, na.rm = T)
      
      EA_sh_RoW <- mean(c$RegType_RoW_EA, na.rm = T)
      DEM_sh_RoW <- mean(c$RegType_RoW_DEM, na.rm = T)
      
      EA_sh_Polity <- mean(c$Politytype_Anocracy, na.rm = T)
      DEM_sh_Polity <- mean(c$Politytype_Democracy, na.rm = T)
      
      EA_sh_FH <- mean(c$status_fh_PF, na.rm = T)
      DEM_sh_FH<- mean(c$status_fh_F, na.rm = T)
      
      EA_sh_AF <- mean(c$AnckarRegtype_MP_Autocracy, na.rm = T)
      DEM_sh_AF<- mean(c$AnckarRegtype_Democracy, na.rm = T)
      
      EA_sh_HTW <- mean(c$HTW_RegType_MP_Autocracy, na.rm = T)
      DEM_sh_HTW<- mean(c$HTW_RegType_Democracy, na.rm = T)
      
      
      EA_sh_AWD <- mean(c$RegType_magaloni_EA, na.rm = T)
      DEM_sh_AWD <- mean(c$RegType_magaloni_DEM, na.rm = T)
      
      
      
      data$neigh_EA_lied <- ifelse(data$year == current_year & data$ccode == current_country, EA_sh_lied, data$neigh_EA_lied )
      data$neigh_DEM_lied <- ifelse(data$year == current_year & data$ccode == current_country, DEM_sh_lied, data$neigh_DEM_lied )
      
      data$neigh_EA_RoW <- ifelse(data$year == current_year & data$ccode == current_country, EA_sh_RoW, data$neigh_EA_RoW )
      data$neigh_DEM_RoW <- ifelse(data$year == current_year & data$ccode == current_country, DEM_sh_RoW, data$neigh_DEM_RoW )
      
      data$neigh_EA_Polity <- ifelse(data$year == current_year & data$ccode == current_country, EA_sh_Polity, data$neigh_EA_Polity )
      data$neigh_DEM_Polity <- ifelse(data$year == current_year & data$ccode == current_country, DEM_sh_Polity, data$neigh_DEM_Polity )
      
      data$neigh_EA_FH <- ifelse(data$year == current_year & data$ccode == current_country, EA_sh_FH, data$neigh_EA_FH )
      data$neigh_DEM_FH <- ifelse(data$year == current_year & data$ccode == current_country, DEM_sh_FH, data$neigh_DEM_FH )
      
      data$neigh_EA_AF <- ifelse(data$year == current_year & data$ccode == current_country, EA_sh_AF, data$neigh_EA_AF )
      data$neigh_DEM_AF <- ifelse(data$year == current_year & data$ccode == current_country, DEM_sh_AF, data$neigh_DEM_AF )
      
      data$neigh_EA_HTW <- ifelse(data$year == current_year & data$ccode == current_country, EA_sh_HTW, data$neigh_EA_HTW )
      data$neigh_DEM_HTW <- ifelse(data$year == current_year & data$ccode == current_country, DEM_sh_HTW, data$neigh_DEM_HTW )
      
      data$neigh_EA_AWD <- ifelse(data$year == current_year & data$ccode == current_country, EA_sh_AWD, data$neigh_EA_AWD )
      data$neigh_DEM_AWD <- ifelse(data$year == current_year & data$ccode == current_country, DEM_sh_AWD, data$neigh_DEM_AWD )
      
      
      
      
      rm(b,c)
      print(paste("Values created for", current_country))
      
      
      
      
      
    }
    
    else{
      print(paste("Current country", current_country, "not existent in", current_year))
      print(paste("Skipping this country for", current_year))
      
      
    }
  }
  
}




summary(data$neigh_DEM_lied)
d2 <- data%>%filter(!is.na(neigh_EA_lied))



save(df, file = "Data/Replication/Kim/Kim_Replication_Ready.rda")








